

!-------------------------------------------------------------------------------------
subroutine evolve_vp_advect(dfu, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, Dzzzz
  use fft_wrap

  implicit none

  real,    intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: dt

  integer                :: iv, ix, iy, iDir
  complex                :: fz(jz)
  complex                :: iVdt
  complex                :: Dkdt(jz)


  !-- advect in z-direction

  iy=jy

  Dkdt(:) =  Dzzzz * kZvec(:)**4 *dt

  do iDir=1,jDir
     do ix=1, jx
        do iv=1,jv
           
           iVdt = iunit * dt * vVec(iv)

           fz(:) = dfu(iv, :, iy, ix, iDir)
           call fft_forward(jz,fz,fz)
           fz = fz/float(jz)
           fz(:) = fz(:) * exp ( -iVdt * kZvec(:) - Dkdt(:) )
           call fft_backward(jz,fz,fz)
           dfu(iv, :, iy, ix, iDir) = fz(:)

        enddo
     enddo
  enddo

end subroutine evolve_vp_advect



!-------------------------------------------------------------------------------------
subroutine evolve_vp_accel(dfu, Ez2D, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, Dvvvv
  use fft_wrap

  implicit none

  real,    intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: Ez2D(jz, jy, jx)
  real,    intent(in)    :: dt

  integer                :: iz, iv, ix, iy, iDir
  complex                :: fv(jv)

  complex                :: iEdt
  complex                :: Dkdt(jv)


  !-- advect in v-direction

  iy = jy

  Dkdt(:) = Dvvvv*kVvec(:)**4*dt 

  do idir=1,jDir
     do ix=1, jx
        do iz=1,jz

           iEdt =  iunit * dt * Ez2D(iz,iy,ix) 

           fv(:) = dfu(:, iz, iy, ix, iDir)
           call fft_forward(jv,fv,fv)
           fv = fv/float(jv)
           fv(:) = fv(:) * exp ( - iEdt * kVvec(:) - Dkdt(:) )
           call fft_backward(jv,fv,fv)
           dfu(:, iz, iy, ix, iDir) = fv(:)

        enddo
     end do
  end do

end subroutine evolve_vp_accel


!-------------------------------------------------------------------------------------
subroutine evolve_vp_accel_lin( dfu, Fbg, Ez2D, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, Dvvvv
  use fft_wrap

  implicit none

  real,    intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: Fbg(jv)
  real,    intent(in)    :: Ez2D(jz, jy, jx)
  real,    intent(in)    :: dt


  complex                :: fv(jv)
  complex                :: expDvvvv(jv), ff0(jv)

  real                   :: dFbg(jv)

  integer                :: ix, iz, iy, iv, iDir


  ix=1
  iy=1
  iDir=1


  !-- derivative of background distribution

  fv(:) = Fbg(:)
  call fft_forward(jv,fv,fv)
  fv = fv/float(jv)

  ff0(:) = iunit*fv(:) / Dvvvv / kVvec(:)**3
  ff0(1) = 0


  fv(:) = iunit*fv(:)*kVvec(:)
  call fft_backward(jv,fv,fv)
  dFbg(:) = fv(:)

  expDvvvv(:) =  exp(-Dvvvv * kVvec(:)**4 *dt)



  !-- advect in v-direction

  if (Dvvvv .eq. 0) then

     do iz=1, jz
        dfu(:, iz, iy, ix, iDir) = dfu(:, iz, iy, ix, iDir) - dFbg(:)*Ez2D(iz,iy,ix)*dt
     end do

  else

     do iz=1,jz

        fv(:) = dfu(:, iz, iy, ix, iDir)
        call fft_forward(jv,fv,fv)
        fv = fv/float(jv)
        fv(:) = fv(:)*expDvvvv + ff0*Ez2D(iz,iy,ix)*(expDvvvv - 1)
        call fft_backward(jv,fv,fv)
        dfu(:, iz, iy, ix, iDir) = fv(:)

     end do

  end if

end subroutine evolve_vp_accel_lin



!-------------------------------------------------------------------------------------
subroutine evolve_vp_efield(dfu, Ez3D)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit
  use fft_wrap

  implicit none

  real,    intent(in)   :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(out)  :: Ez3D(jz, jy, jx)

  complex               :: fz(jz)
  integer               :: ix, iz, iy, iv, idir



  !-- compute charge density
 
  ix=1
  iy=1
 
  do iz=1,jz
     fz(iz) = sum(dfu(:,iz,iy,ix,:))*dv
  end do
 

  !-- find electric field --

  call fft_forward(jz,fz,fz)
  fz = fz/float(jz)
  fz = -iunit*fz/kZvec
  fz(1) = 0 
  call fft_backward(jz,fz,fz)

  Ez3D(:,1,1) = fz


end subroutine evolve_vp_efield


!-------------------------------------------------------------------------------------
subroutine evolve_vp_filter(dfu)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: Vmax, Vcutoff


  implicit none

  real, intent(inout)   :: dfu(jv, jz, jy, jx, jDir)

  integer               :: iv
  logical, save         :: first_call = .true.

  real, save            :: filter(jv)


  if (Vcutoff .eq. 0) return


  !-- generate filter --

  if (first_call) then

     if (Vcutoff .gt. 0) then
        filter(:) = vvec(:) / (Vcutoff*Vmax)
        filter(:) = exp( - filter(:)**(16) )
     else
        filter(:) = 1 - (0.5 + 0.5*tanh( (abs(vvec(:)/Vmax) + Vcutoff)*16))**4
     end if

     first_call = .false.

  end if



  !-- apply filter --
  
  do iv=1,jv

    dfu(iv,:,:,:,:) = dfu(iv,:,:,:,:)*filter(iv)

  end do


end subroutine evolve_vp_filter


!-------------------------------------------------------------------------------------



